function [rkstat, pval_Wald, pval_F]=rank_eig(R, F)
%{
    This function implements the rank test (eigenvalue approach)
    
    INPUT 
    R       : n by T, a matrix of asset returns
    F       : k by T, a matrix of risk factors
    
    OUTPUT  
    rkstat  : rank test statistic


    NOTATION: |\theta*B-A|=0, see Kleibergen (2009)
%}

    [n, T]=size(R);                       % get n, T
    k=size(F,1);                          % get k
    F_bar=F-mean(F,2)*ones(1,T);          % derive F_bar
    B=inv(F_bar*F_bar'/T);
     
    ln=ones(n-1,1);                       % a vector of ones
    Rs=R(2:n,:)-kron(ln,mean(R(:,:)));
    %Rs=R(1:(n-1),:)-kron(ln,R(n,:));      % delete the n_th asset    
    F_exp=[F;ones(1,T)];
    B_temp=Rs*F_exp'/(F_exp*F_exp');
    Sigma_tt=1/T*(Rs-B_temp*F_exp)*(Rs-B_temp*F_exp)';
    Bhat=B_temp(:,1:end-1);
    A=Bhat'*inv(Sigma_tt)*Bhat;           % matrix A
    
    % rank test
    [V,D]=eig(A,B);                       % V:eigenvectors; D:eigenvalues  
    
    V
    [mval,mind]=min(diag(D));
    mind
    V(:,mind)' * Bhat'
    
    
    rkstat=T*min(diag(D));
    pval_Wald=1-chi2cdf(rkstat,n-k);
    cRank=k-1;
    dg1=(n-1-cRank)*(k-cRank);	
    dg2=(T-k-1)-dg1+1;
    stat_F=rkstat/T*dg2/dg1;
    pval_F=1-fcdf(stat_F,dg1,dg2);
end 